classdef SLM < handle
    properties
        deviceModel
        deviceWidth
        deviceHeight
        pixelSize
        monitorPos
        wavelength
        wavenumber
        focallength
        
        correction
        coordinate
        transform
        
        deviceFigure
        deviceAxes
        
        totalIntensity
        incIntensity
        targetIntensity
        targetNonzeroIndex
        maskMat
    end
    
    methods
        function obj = SLM(varargin)
            p = inputParser;
            p.addParameter('waveLength', 780, @(x) x>=750 && x<=850); % nm
            p.addParameter('focalLength', 0.3, @(x) x>0); %m
            p.parse(varargin{:});
            
            obj.deviceModel = 'X13138';
            obj.deviceWidth = 1280;
            obj.deviceHeight = 1024;
            obj.pixelSize = 12.5e-6;% m
            
            obj.findDevice();

            obj.wavelength = p.Results.waveLength;% nm
            obj.wavenumber = 2.0*pi/(obj.wavelength*1e-9); % m^{-1}
            obj.focallength = p.Results.focalLength; %m

            obj.coordinate = obj.defineInputPlane();
            obj.transform = obj.defineFourierPlane();
            
            obj.newPlot();
            obj.loadCorrectionFile(obj.wavelength);
            obj.constPhase();
        end
        
        function delete(obj)
            L = log4m.getLogger;
            try
                close(obj.deviceFigure);
                L.info('EquipmentControl.Hamamatsu.SLM::delete', 'Current Figure is closed.');
            catch
                L.info('EquipmentControl.Hamamatsu.SLM::delete', 'Current Figure not found.');
            end
        end
        
        function findDevice(obj)
            L = log4m.getLogger;
            mp = get(0, 'MonitorPositions');
            nMonitor = size(mp, 1);
            if nMonitor == 1
                L.info('EquipmentControl.Hamamatsu.SLM::findDevice', 'Only one monitor is detected.');
                obj.monitorPos = [1, 1, obj.deviceWidth, obj.deviceHeight];
            else
                idx = find(mp(:, 3) == obj.deviceWidth & mp(:, 4) == obj.deviceHeight);
                if isempty(idx)
                    L.info('EquipmentControl.Hamamatsu.SLM::findDevice', sprintf('None of [%d] monitors is fitted to the device.', nMonitor));
                    obj.monitorPos = [1, 1, obj.deviceWidth, obj.deviceHeight];
                elseif length(idx) == 1
                    L.info('EquipmentControl.Hamamatsu.SLM::findDevice', sprintf('One monitor is detected: [left, bottom] = [%d, %d].', mp(idx, 1:2)) );
                    obj.monitorPos = mp(idx, :);
                else
                    L.info('EquipmentControl.Hamamatsu.SLM::findDevice', sprintf('%d monitors are detected. Using the first one: [left, bottom] = [%d, %d].', length(idx), mp(idx, 1:2)) );
                    obj.monitorPos = mp(idx(1), :);
                end
            end
        end
        
        function coord = defineInputPlane(obj)
            coord.xList = (1:obj.deviceWidth)-obj.deviceWidth/2;
            coord.yList = (1:obj.deviceHeight)-obj.deviceHeight/2;
            [coord.xMat, coord.yMat]=meshgrid(coord.xList, coord.yList);
            coord.zMat=coord.xMat+1j*coord.yMat;
            coord.rMat=abs(coord.zMat);
            coord.pMat=angle(coord.zMat);
        end
        
        function transform = defineFourierPlane(obj)
            transform.kxList = (-obj.deviceWidth/2:obj.deviceWidth/2-1)*(2*pi/obj.pixelSize/obj.deviceWidth);
            transform.kyList = (-obj.deviceHeight/2:obj.deviceHeight/2-1)*(2*pi/obj.pixelSize/obj.deviceHeight);
            [transform.kxMat, transform.kyMat]=meshgrid(transform.kxList, transform.kyList);
            
            transform.xPixelSize = pi/obj.pixelSize/obj.deviceWidth*obj.focallength/obj.wavenumber;
            transform.yPixelSize = pi/obj.pixelSize/obj.deviceHeight*obj.focallength/obj.wavenumber;
            transform.xList = transform.kxList*obj.focallength/obj.wavenumber;
            transform.yList = transform.kyList*obj.focallength/obj.wavenumber;
            [transform.xMat, transform.yMat]=meshgrid(transform.xList, transform.yList);
            transform.zMat=transform.xMat+1j*transform.yMat;
            transform.rMat=abs(transform.zMat);
            transform.pMat=angle(transform.zMat);
        end
        
        function loadCorrectionFile(obj, varargin)
            L = log4m.getLogger;
            p=inputParser;
            p.addOptional('wavelength', obj.wavelength, @(x) ismember(x, 750:10:850)); % nm
            p.parse(varargin{:});
            
            d = Constant.PathDelimiter;
            filename = sprintf('%s%c+EquipmentControl%c+Hamamatsu%c@SLM%cdeformation_correction_pattern%cCAL_LSH0801993_%dnm.bmp', ...
                            Constant.AutoMeasurementCorePath, d, d, d, d, d, p.Results.wavelength);
            filedata = imread(filename);
            
            c.Image = [repmat(filedata(:, 1), 1, 4) filedata repmat(filedata(:, end), 1, 4)];
            c.Phase = mat2gray(c.Image, [0 255])*2*pi-pi;
            c.Phase1 = unwrap(unwrap(c.Phase), [], 2);
            obj.correction = c;
            
            L.info('EquipmentControl.Hamamatsu.SLM::loadCorrectionFile', sprintf('Deformation Correction Pattern of %d nm loaded.', p.Results.wavelength));
        end
        
        function newPlot(obj)
            obj.deviceFigure = figure('MenuBar', 'none', 'ToolBar', 'none', 'Position', obj.monitorPos);
            obj.deviceAxes = axes(obj.deviceFigure, 'Units', 'normalized', 'Position', [0 0 1 1]);
        end
        
        function dispImage(obj, phaseMat, varargin)
            if ~isvalid(obj.deviceFigure) || ~isvalid(obj.deviceAxes)
                obj.newPlot();
            end
            
            imageData = obj.phase2image(phaseMat);            
            imshow(imageData, 'Parent', obj.deviceAxes, 'colormap', colormap('gray'));            
        end
        
        function imageData = phase2image(~, phaseMat)
            imageData = uint8(mat2gray(phaseMat, [-pi pi])*256);
        end
        
        %%        
        function phaseMat = dcp(obj, wavelength) % deformation correction pattern
            if nargin == 2 % overwrite default wavelength
                obj.wavelength = wavelength;
                obj.loadCorrectionFile(wavelength);
            end
            phaseMat = obj.correction.Phase;
        end
        
        function phaseMat = shift(obj, distX, distY, distZ)
            kx = distX * obj.wavenumber/obj.focallength;
            ky = distY * obj.wavenumber/obj.focallength;
            kz = distZ * obj.wavenumber/obj.focallength;
            complexZmat = exp(1j.* kx.*obj.pixelSize.*obj.coordinate.xMat ...
                            + 1j.* ky.*obj.pixelSize.*obj.coordinate.yMat ...
                            + 1j.* kz.*(obj.pixelSize.*obj.coordinate.rMat).^2)./(2.0*obj.focallength);
            phaseMat = angle(complexZmat);
        end
        
        function phaseMat = constPhase(obj, varargin)
            p=inputParser;
            p.addOptional('value', 0.0, @isnumeric);
            p.addParameter('isDCP', true, @islogical);
            p.addParameter('isDisp', true, @islogical);
            p.parse(varargin{:});
            
            phaseMat = p.Results.value.*ones(size(obj.coordinate.xMat));
            
            if p.Results.isDCP
                phaseMat = phaseMat + obj.correction.Phase;
            end
            
            if p.Results.isDisp
                obj.dispImage(phaseMat);
            end
        end

        function phaseMat = dot(obj, varargin)
            p=inputParser;
            p.addOptional('distX', 0.0, @isnumeric);
            p.addOptional('distY', 0.0, @isnumeric);
            p.addOptional('distZ', 0.0, @isnumeric);
            p.addParameter('isDCP', true, @islogical);
            p.addParameter('isDisp', true, @islogical);
            p.parse(varargin{:});
            
            phaseMatShift = obj.shift(p.Results.distX, p.Results.distY, p.Results.distZ);
            phaseMat = phaseMatShift;
            
            if p.Results.isDCP
                phaseMat = phaseMat + obj.correction.Phase;
            end
            
            if p.Results.isDisp
                obj.dispImage(phaseMat);
            end
        end
        
        function phaseMat = line(obj, len, w0, varargin)
            %Example: slm.line(0.0005, 20e-6, 0.00, 0.001, 'dir', [1 1]);
            p=inputParser;
            p.addOptional('distX', 0.0, @isnumeric);
            p.addOptional('distY', 0.0, @isnumeric);
            p.addOptional('distZ', 0.0, @isnumeric);            
            p.addParameter('dir', [1, 0], @(x) isnumeric(x) && length(x)==2);
            p.addParameter('isDCP', true, @islogical);
            p.addParameter('isDisp', true, @islogical);
            p.parse(varargin{:});
            
            cL = obj.wavenumber^2 * w0*len /(8.0*obj.focallength^2)*obj.pixelSize^2;% cL = k^2 * w0 * L /8 /f^2 = k/(2F) Alxe Gaunt Table 5.1
            normX = p.Results.dir(1)/norm(p.Results.dir);
            normY = p.Results.dir(2)/norm(p.Results.dir);
            coordMat2 = (normX*obj.coordinate.xMat + normY*obj.coordinate.yMat).^2;
            complexZmat = exp(-1j.* cL.* coordMat2);
            
            phaseMatShift = obj.shift(p.Results.distX, p.Results.distY, p.Results.distZ);
            phaseMat = angle(complexZmat) + phaseMatShift;
            
            if p.Results.isDCP
                phaseMat = phaseMat + obj.correction.Phase;
            end
            
            if p.Results.isDisp
                obj.dispImage(phaseMat);
            end
        end
        
        function phaseMat = vortex(obj, l, varargin)
            p=inputParser;
            p.addOptional('z0', 0.0, @isnumeric);
            p.addOptional('distX', 0.0, @isnumeric);
            p.addOptional('distY', 0.0, @isnumeric);
            p.addOptional('distZ', 0.0, @isnumeric);
            p.addParameter('isDCP', true, @islogical);
            p.addParameter('isDisp', true, @islogical);
            p.parse(varargin{:});

            phaseMatShift = obj.shift(p.Results.distX, p.Results.distY, p.Results.distZ);
            phaseMat = angle((obj.coordinate.zMat-p.Results.z0).^l) + phaseMatShift;
            
            if p.Results.isDCP
                phaseMat = phaseMat + obj.correction.Phase;
            end
            
            if p.Results.isDisp
                obj.dispImage(phaseMat);
            end
        end
        
        function phaseMat = axicon(obj, r, varargin)
            p=inputParser;
            p.addOptional('distX', 0.0, @isnumeric);
            p.addOptional('distY', 0.0, @isnumeric);
            p.addOptional('distZ', 0.0, @isnumeric);
            p.addParameter('isDCP', true, @islogical);
            p.addParameter('isDisp', true, @islogical);
            p.parse(varargin{:});
            
            
            complexZmat = exp(-1j.* obj.wavenumber.*r./obj.focallength .* obj.pixelSize.*obj.coordinate.rMat);
            
            phaseMatShift = obj.shift(p.Results.distX, p.Results.distY, p.Results.distZ);
            phaseMat = angle(complexZmat) + phaseMatShift;
            
            if p.Results.isDCP
                phaseMat = phaseMat + obj.correction.Phase;
            end
            
            if p.Results.isDisp
                obj.dispImage(phaseMat);
            end
        end

        %%
        function setInputGaussian(obj, beamwidth)
            obj.incIntensity = exp(- (obj.pixelSize.*obj.coordinate.rMat).^2/2/beamwidth^2);
            obj.totalIntensity = sqrt(sum(abs(obj.incIntensity(:)).^2));
        end
        
        function setTargetIntensity(obj, mat)
            normfact = sqrt(sum(abs(mat(:)).^2));
            obj.targetIntensity = abs(mat)./normfact.*obj.totalIntensity;
            obj.targetNonzeroIndex = find(obj.targetIntensity>1./256);
        end
        
        function setTargetFlatTop(obj, shape, param)
            switch lower(shape)
                case 'square'
                    try 
                        width = 0.5*param.width;
                        height = 0.5*param.height;
                        order = param.order;
                    catch
                        error('Parameter is incorrect.');
                    end
                    funcX = 1./( (obj.transform.xMat/width).^order+1);
                    funcY = 1./( (obj.transform.yMat/height).^order+1);
                    mat = funcX.*funcY; normfact = sqrt(sum(abs(mat(:)).^2));
                    obj.targetIntensity = abs(mat)./normfact.*obj.totalIntensity;
                    obj.targetNonzeroIndex = find(obj.targetIntensity>1./256);
                case 'circle'
                    try 
                        radius = param.radius;
                        order = param.order;
                    catch
                        error('Parameter is incorrect.');
                    end
                    mat = 1./( (obj.pixelSize*obj.coordinate.rMat/radius).^order+1 );
                    normfact = sqrt(sum(abs(mat(:)).^2));
                    obj.targetIntensity = abs(mat)./normfact.*obj.totalIntensity;
                    obj.targetNonzeroIndex = find(obj.targetIntensity>1./256);
                otherwise
            end
        end
        
        function showTarget(obj)
            figure;
            imagesc(obj.transform.xList, obj.transform.yList, obj.targetIntensity);
        end
        
        function setMaskRegion(obj, xRange, yRange)
            obj.maskMat = zeros(obj.deviceHeight, obj.deviceWidth);
            obj.maskMat(obj.deviceHeight/2-yRange:obj.deviceHeight/2+yRange, obj.deviceWidth/2-xRange:obj.deviceWidth/2+xRange)=1;
        end
        
        function setMaskRegionShape(obj, shape, param)
            switch lower(shape)
                case 'square'
                    try 
                        width = 0.5*param.mask_width;
                        height = 0.5*param.mask_height;
                    catch
                        error('Parameter is incorrect.');
                    end
                    obj.maskMat = zeros(obj.deviceHeight, obj.deviceWidth);
                    obj.maskMat( abs(obj.transform.xMat)<width & abs(obj.transform.yMat)<height ) = 1;
                case 'circle'
                    try
                        radius = param.mask_radius;
                    catch
                        error('Parameter is incorrect.');
                    end
                    obj.maskMat = zeros(obj.deviceHeight, obj.deviceWidth);
                    obj.maskMat( obj.transform.rMat<radius ) = 1;
                otherwise
            end
        end
        function showMask(obj)
            figure;
            imagesc(obj.transform.xList, obj.transform.yList, obj.maskMat);
        end
        
        function [phaseMatFocalPlane, intensityMatFocalPlane] = dispTransform(obj, phaseMat, varargin)
            L = log4m.getLogger;
            p=inputParser;
            p.addParameter('figure_handle', gcf, @(x) isa(x,  'matlab.ui.Figure'));
            p.addParameter('xlim', [-3, 3], @(x) isnumeric(x) && length(x)==2 && x(2) > x(1));
            p.parse(varargin{:});
            
            ax1 = subplot(2, 3, 1, 'Parent', p.Results.figure_handle); 
            imagesc(ax1, obj.incIntensity); colormap(ax1, 'gray'); title(ax1, 'Input Intensity'); colorbar(ax1);
            xlim(ax1, [0 obj.deviceWidth]); ylim(ax1, [0 obj.deviceHeight]); 
            xlabel(ax1, 'Pixel Number');  ylabel(ax1, 'Pixel Number');
            
            ax4 = subplot(2, 3, 4, 'Parent', p.Results.figure_handle);
            imagesc(ax4, phaseMat); colormap(ax4, 'gray'); title(ax4, 'SLM Phase'); colorbar(ax4);
            xlim(ax4, [0 obj.deviceWidth]); ylim(ax4, [0 obj.deviceHeight]);
            xlabel(ax4, 'Pixel Number'); ylabel(ax4, 'Pixel Number');
            
            [phaseMatFocalPlane, intensityMatFocalPlane] = obj.forwardPropagate( obj.incIntensity.*exp(1j.*phaseMat));            
            [eff, ~] = obj.evaluation(intensityMatFocalPlane);
            
            ax2 = subplot(2, 3, 2, 'Parent', p.Results.figure_handle);
            imagesc(ax2, obj.transform.xList*1e3, obj.transform.yList*1e3, intensityMatFocalPlane); colormap(ax2, 'jet');
            aspectratio = pbaspect(ax2);
            xlim(ax2, p.Results.xlim); ylim(ax2, p.Results.xlim*aspectratio(2));
            xlabel(ax2, 'Position (mm)'); ylabel(ax2, 'Position (mm)');
            title(ax2, 'Transformed Intensity'); colorbar(ax2);

            ax3 = subplot(2, 3, 3, 'Parent', p.Results.figure_handle);
            plot(ax3, obj.transform.yList*1e3, intensityMatFocalPlane(:, obj.deviceWidth/2), 'r-'); hold(ax3, 'on');
            try
                plot(ax3, obj.transform.yList*1e3, obj.targetIntensity(:, obj.deviceWidth/2)*eff, 'b-');
            catch
                L.info('EquipmentControl.Hamamatsu.SLM::dispTransform', 'Target intensity not found.');
            end
             xlim(ax3, p.Results.xlim); title(ax3, 'Vertical Cut');
             xlabel(ax3, 'Position (mm)'); ylabel(ax3, 'Intenstity');
             view(ax3, [90 90]); grid(ax3, 'on'); hold(ax3, 'off');

            ax5 = subplot(2, 3, 5, 'Parent', p.Results.figure_handle);
            plot(ax5, obj.transform.xList*1e3, intensityMatFocalPlane(obj.deviceHeight/2, :), 'r-'); hold(ax5, 'on');
            try
                plot(ax5, obj.transform.xList*1e3, obj.targetIntensity(obj.deviceHeight/2, :)*eff, 'b-');
            catch
                L.info('EquipmentControl.Hamamatsu.SLM::dispTransform', 'Target intensity not found.');
            end
             xlim(ax5, p.Results.xlim); title(ax5, 'Horizontal Cut');
             xlabel(ax5, 'Position (mm)'); ylabel(ax5, 'Intenstity');
             grid(ax5, 'on'); hold(ax5, 'off');
            
            ax6 = subplot(2, 3, 6, 'Parent', p.Results.figure_handle);
            imagesc(ax6, obj.transform.xList*1e3, obj.transform.yList*1e3, phaseMatFocalPlane); colormap(ax6, 'gray');
            aspectratio = pbaspect(ax6);
            xlim(ax6, p.Results.xlim); ylim(ax6, p.Results.xlim*aspectratio(2));
            xlabel(ax6, 'Position (mm)'); ylabel(ax6, 'Position (mm)');
            title(ax6, 'Transformed Phase'); colorbar(ax6);
        end
        
        function [phaseMat, ampMat] = forwardPropagate(obj, zMatIn)
            zMat = fftshift(fft2(zMatIn));
            phaseMat = angle(zMat);
            ampMat1 = abs(zMat);
            normFactor = sqrt(sum(ampMat1(:).^2));
            ampMat = abs(zMat)./normFactor.*obj.totalIntensity;
        end
        function [phaseMat, ampMat] = backwardPropagate(obj, zMatOut)
            zMat = ifft2(ifftshift(zMatOut));
            phaseMat = angle(zMat);
            ampMat1 = abs(zMat);
            normFactor = sqrt(sum(ampMat1(:).^2));
            ampMat = abs(zMat)./normFactor.*obj.totalIntensity;
        end
        
        function mix = amplitudeMixing(obj, ampMat, m, delta)
            mat1 = m*obj.maskMat.*sqrt(obj.targetIntensity.^2 + delta.^2);
            mat2 = (1-m)* (1-obj.maskMat).*ampMat;
            mix = mat1./max(mat1(:)) + mat2./max(mat2(:));
        end
        
        function [efficent, rms] = evaluation(obj, intensityMat_FP)
            efficent = sqrt(sum(intensityMat_FP(obj.targetNonzeroIndex).^2))/obj.totalIntensity;
            rms = sqrt(sum((intensityMat_FP(obj.targetNonzeroIndex)/efficent - obj.targetIntensity(obj.targetNonzeroIndex)).^2))/obj.totalIntensity;
        end
        
        function [phaseMat_SLM, intensityMat_FP, phaseMat_FP, efficient, rms] = MRAF(obj, iter, m, delta, varargin)
            p = inputParser;
            p.addParameter('isPlotProgress', false, @islogical);
            p.addParameter('isPlotConvergence', true, @islogical);
            p.parse(varargin{:});
            
            clear textprogressbar
            
            if p.Results.isPlotProgress
                f = figure; 
                mm(iter) = struct('cdata',[],'colormap',[]);
                
            end
            
            if p.Results.isPlotConvergence
                f1 = figure; ax = axes('Parent', f1);
                effList = nan(1, iter);
                rmsList = nan(1, iter);
            end
            
            zMatSLMplane = obj.incIntensity .* exp(1j.*obj.backwardPropagate(obj.targetIntensity));
            textprogressbar('Iteration progress: ');
            for k =1:iter
                textprogressbar(k/iter*100);
                
                [phaseMat_FP, ampMat_FP] = obj.forwardPropagate(zMatSLMplane);
                fourierPlaneIntensity = obj.amplitudeMixing(ampMat_FP, m, delta);
                zMatFourierPlane = fourierPlaneIntensity .* exp(1j.*phaseMat_FP);
                phaseMat_SLM = obj.backwardPropagate(zMatFourierPlane);
                zMatSLMplane =  obj.incIntensity .* exp(1j.*phaseMat_SLM);
                
                if p.Results.isPlotConvergence
                    [effList(k), rmsList(k)] = obj.evaluation(ampMat_FP);
                    yyaxis(ax, 'left'); semilogy(ax, 1:iter, effList, 'o-'); ylabel(ax, 'Efficiency'); 
                    yyaxis(ax, 'right'); semilogy(ax, 1:iter, rmsList, 'd-'); ylabel(ax, 'Relative deviation');
                    xlabel(ax, 'iteration number'); xlim(ax, [1, iter]); grid(ax, 'on'); legend(ax, {'Efficiency', 'Deviation'}, 'Location', 'best');
                    drawnow; 
                end
                    
                if p.Results.isPlotProgress
                    obj.dispTransform(phaseMat_SLM, 'figure_handle', f); drawnow;                    
                    mm(k) = getframe(f);
                end
            end
            textprogressbar('\t Done.');
            
            figure;
            [phaseMat_FP, intensityMat_FP] = obj.dispTransform(phaseMat_SLM);
            [efficient, rms] = obj.evaluation(intensityMat_FP);
        end
    end
end

